function [betas, irtc, ptc, se, weights] = sdm_CalcBetas(hfile, tcd, dim)
% SDM::CalcBetas  - perform GLM calculcation
%
% FORMAT:       [betas, irtc, ptc, se] = sdm.CalcBetas(tc [, tdim])
%
% Input fields:
%
%       tc          time course data (numeric data or FMR/VTC/MTC)
%       tdim        required for numeric tc data, temporal dimension
%                   alternatively, a struct can be given with options
%        .robust    perform robust regression instead of OLS
%        .tdim      temporal dimension (otherwise auto-detect)
%        .tssm      either SSM or TSM object with mappings
%
% Output fields:
%
%       betas       betas, in time course data dimension
%       irtc        inverse design matrix
%       ptc         predicted tc
%       se          standard error
%       weights     robust-regression weights
%
% Note: if no constant term is in the matrix, it will be added as the
%       last column
%       also, for robust regression, the standard error, if requested
%       is adapted to the nominal d.f.

% Version:  v0.8a
% Build:    9102000
% Date:     Oct-20 2009, 12:46 AM CEST
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'sdm')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end

% get design matrix right
bc = bvqxfile_getcont(hfile.L);
sdm = bc.SDMMatrix;

% further arguments
robust = false;
tssm = [];
if nargin > 2 && ...
    isstruct(dim) && ...
    numel(dim) == 1
    if isfield(dim, 'robust') && ...
        islogical(dim.robust) && ...
        numel(dim.robust) == 1
        robust = dim.robust;
    end
    if isfield(dim, 'tssm') && ...
        numel(dim.tssm) == 1 && ...
       (isBVQXfile(dim.tssm, 'ssm') || ...
        isBVQXfile(dim.tssm, 'tsm')) && ...
        numel(tcd) == 1 && ...
        isBVQXfile(tcd, 'mtc')
        tssm = dim.tssm;
    end
    if isfield(dim, 'tdim') && ...
        isa(dim.tdim, 'double') && ...
        numel(dim.tdim) == 1 && ...
        any((1:5) == dim.tdim)
        dim = dim.tdim;
    end
end
if nargin < 3 || ...
   ~isa(dim, 'double') || ...
    numel(dim) ~= 1 || ...
    isinf(dim) || ...
    isnan(dim) || ...
    fix(dim) ~= dim || ...
    dim < 1 || ...
    dim > 5
    rdim = [];
else
    rdim = dim;
end

% add all-1 column for confound
mrtc = mean(sdm);
vrtc = var(sdm);
if ~any(mrtc ~= 0 & vrtc == 0) && ...
    size(sdm, 2) < size(sdm, 1)
    sdm(:, end + 1) = 1;
end

% number of rows and betas
numrows = size(sdm, 1);
numbets = size(sdm, 2);

% get betas that can be ~= 0
usebets = (sum(abs(sdm)) ~= 0);

% check second argument -> BVQXfile object
if numel(tcd) == 1 && ...
    isBVQXfile(tcd, true)

    % get struct version
    sstc = bvqxfile_getscont(tcd.L);
    stcd = sstc.C;

    % get filetype
    try
        ftype = lower(sstc.S.Extensions{1});
    catch
        error( ...
            'BVQXfile:InternalError', ...
            'Invalid BVQXfile object given as time course data.' ...
        );
    end

    % switch over filetype
    switch lower(ftype)

        % FMR
        case {'fmr'}
            
            % rdim MUST be 3
            rdim = 3;

            % STCs loaded?
            if isempty(stcd.Slice) || ...
               ~isstruct(stcd.Slice) || ...
               ~isfield(stcd.Slice, 'STCData')
                try
                    fmr_LoadSTCs(tcd);
                    stcd = bvqxfile_getcont(tcd.L);
                catch
                    error( ...
                        'BVQXfile:InternalError', ...
                        'Error loading STCs for FMR %s.', ...
                        sstc.F ...
                    );
                end
            end

            % check MTC
            nvol = numel(stcd.Slice(1).STCData(1, 1, :, 1));
            if nvol ~= numrows
                error( ...
                    'BVQXfile:InternalError', ...
                    'Invalid number of time points in FMR/STC.' ...
                );
            end

            % get STC data
            switch (stcd.DataStorageFormat)
                
                % old FMR format (with multiple STCs)
                case {1}
                    nslc = stcd.NrOfSlices;
                    xres = stcd.ResolutionX;
                    yres = stcd.ResolutionY;
                    tcd = zeros(xres, yres, nvolc, nslc);
                    for sc = 1:nslc
                        tcd(:, :, :, sc) = stcd.Slice(sc).STCData(:, :, :);
                    end
                    
                % new FMR format (one STC file)
                case {2}
                    tcd = double(stcd.Slice.STCData(:, :, :, :));
                    
                % unknown FMR format (?)
                otherwise
                    error( ...
                        'BVQXfile:InvalidObject', ...
                        'Unsupported FMR DataStorageFormat.' ...
                    );
            end

        % MTC
        case {'mtc'}

            % rdim MUST be 1
            rdim = 1;
            
            % check MTC
            if size(stcd.MTCData, 1) ~= numrows
                error( ...
                    'BVQXfile:InternalError', ...
                    'Invalid number of time points in MTC.' ...
                );
            end

            % get MTC data
            tcd = double(stcd.MTCData(:, :));
            
            % apply SSM/TSM ?
            if ~isempty(tssm)
            end

        % VTC
        case {'vtc'}

            % rdim MUST be 1
            rdim = 1;
            
            % check VTC
            if size(stcd.VTCData, 1) ~= numrows
                error( ...
                    'BVQXfile:InternalError', ...
                    'Invalid number of volumes in VTC.' ...
                );
            end

            % get VTC data
            tcd = double(stcd.VTCData(:, :, :, :));

        % otherwise
        otherwise

            error( ...
                'BVQXfile:TypeNotSupported', ...
                'BVQXfile object of type %s not supported.', ...
                ftype ...
            );
    end

% for other numerics
elseif isnumeric(tcd) && ...
   ~isa(tcd, 'double')

    % create double
    tcd = double(tcd);

end

% check dim
if isempty(rdim)
    rdim = findfirst(size(tcd) == size(sdm, 1));
    if isempty(rdim)
        error( ...
            'BVQXfile:InternalError', ...
            'Invalid data matrix size vs. given data.' ...
        );
    end
elseif rdim < 1 || ...
    rdim > length(size(tcd)) || ...
    size(tcd, rdim) ~= numrows
    error( ...
        'BVQXfile:InternalError', ...
        'Invalid data matrix vs. selected dimension.' ...
    );
end
if rdim > 1
    neword = [rdim, 1:length(size(tcd))];
    newodr = find(neword == rdim);
    neword(newodr(2)) = [];
    tcd = permute(tcd, neword);
    [newsrt, oldord] = sort(neword);
else
    oldord = [];
end

% reshaping data to comply
tcds = size(tcd);
numvox = prod(tcds(2:end));
tcd = reshape(tcd, [tcds(1), numvox]);
tcds(1) = [];
tcdrs = tcds;
if length(tcdrs) < 2
    tcdrs(2) = 1;
end

% robust
if robust
    switch (nargout)
        case {1, 2}
            betas = fitrobustbisquare_img(sdm(:, usebets), tcd');
        case {3, 4}
            [betas, weights] = fitrobustbisquare_img(sdm(:, usebets), tcd');
            weights = weights';
            ptc = (sdm(:, usebets) * betas') .* weights + (tcd .* (1 - weights));
            if nargout > 3
                se = squeeze(sqrt((1 / (numrows - numbets)) .* (sum(weights) - 1)) .* ...
                    std(tcd - ptc, 0));
                se(isinf(se) | isnan(se)) = Inf;
            end
    end
    
% or OLS
else
    switch (nargout)
        case {1}
            betas = calcbetas(sdm(:, usebets), tcd, 1);
        case {2}
            [betas, irtc] = calcbetas(sdm(:, usebets), tcd, 1);
        case {3}
            [betas, irtc, ptc] = calcbetas(sdm(:, usebets), tcd, 1);
        case {4}
            [betas, irtc, ptc, se] = calcbetas(sdm(:, usebets), tcd, 1);
    end
    if nargout > 4
        weights = ones(size(tcd));
    end
end

% not all betas?
if sum(usebets) ~= numbets
    usebetd = ones(size(usebets));
    usebetd(usebets) = 1:sum(usebets);
    betas = betas(:, usebetd);
    betas(:, ~usebets) = 0;
    
    % also build standard irtc ?
    if nargout > 1
        irtc = pinv(sdm' * sdm);
        irtc(~usebets, :) = 0;
        irtc(:, ~usebets) = 0;
    end
end

% reshape betas
betas = reshape(betas, [tcds, numbets]);

% reshape ptc if necessary
if nargout > 2
    ptc = reshape(ptc, [numrows, tcds]);
    if ~isempty(oldord)
        if nargout > 4
            weights = reshape(weights, size(ptc));
            weights = permute(weights, oldord);
        end
        ptc = permute(ptc, oldord);
    end
    if nargout > 3
        se = reshape(se, tcdrs);
    end
end
